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<^ ' Parameters of the nonorthogonal tight-binding model for hydrocarbons are derived 

based on a criterion of the best agreement between the calculated and experimental 

values of bond lengths and binding energies for different molecules C^H^. The results 
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+JDbtained can be used, e. g., to study the kinetics of hydrogen absorption by carbon 



inianostructures, to simulate the dynamics of hydrocarbon clusters like cubane CgHg, etc. 
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1. Introduction 

In recent years, there was a growing interest in interaction between hydrogen (both 
atomic and molecular) and various carbon structures. This is in part due to the promising 
results of the first experiments on hydrogen absorption by carbon nanotubes P, 121 E] and 
a corresponding perspective to use hydrocarbons in hydrogen energetics, see reviews 
[H [3] . A solid based on clusters CgHg (cubanes) [6] is one more example of hydrocarbon 
systems important from both fundamental and practical viewpoints. Cubane and its 
derivatives are considered as candidates to a new type of fuel, they can probably be used 
in pharmaceutics, liquid crystals, etc. [7]. 

In the case that experimental data remain incomplete and controversial, theoretical 
calculations become of particular importance. For simulations on hydrocarbons, either 
classical interatomic potentials, e. g., Tersoff-Brenner potential P, E], or so called first 
principle methods, e. g., the density functional theory (DFT) pHl [TT| . are commonly 



employed. One should keep in mind, however, that making use of ab initio approaches 
do not guarantee, by itself, the reliable results since those approaches are sometimes 
extremely sensitive to the type of exchange-correlation potential, the set of the basis 
functions, the cutoff energy, etc. (for example, different modifications of DFT disagree on 
the most stable C20 isomer, see references in |I2])- Besides, the first principle calculations 
are computer-time consuming and put severe restrictions on the size of the system and/or 
the real-time interval during which the system evolution can be followed (usually ~ 1 — 10 
ps). 

In such a situation, tight-binding models are a reasonable compromise between more 
rigorous ab initio and oversimplified empirical calculations. Contrary to classical poten- 
tials, those models exphcitly account for the contribution of electron subsystem to the 
total energy. Although the tight-binding method is not as strict as ab initio approaches, 



it is competitive with them in the accuracy and not computer-resource-intensive, thus fa- 
cihtating the simulation of relatively large systems and/or dynamic processes at real-time 
scales up to ~ 1 fis. For carbon systems, various tight-binding models were suggested 
in Refs. |I31, HH US, dSl They were successfully used by many authors to simu- 



late both carbon clusters and bulks. In particular, previously we made use of a model 
to study the thermal stability of fuUerenes C20 and Ceo, one-dimensional chains and 
two-dimensional complexes of fullerenes C20, etc. |[T5l [IS, 1221 1211 1221 123, 1211 125 • 



Different tight-binding models for hydrocarbons were suggested in Refs. [221 1271 [2H 



2Sl EDI EH E21- Each of those models has its own advantages and drawbacks, but all of 



them rather accurately describe a broad range of CnHm molecules and macroscopic C-H 
systems (for example, hydrogen on the diamond surface). Among the papers mentioned 
above, of particular interest is Ref. [32] whose authors generalized the tight-binding 



model to C-H-0 systems containing oxygen along with carbon and hydrogen. The model 
|E2] is nonorthogonal (this being important for description of systems with different coor- 
dination numbers [33]) and free of Hubbard-like terms, thus avoiding the self-consistent 
calculation of occupation numbers and reducing the computation time. This model gives 
the geometrical and energetical characteristics of small CnHm and CnHniOi molecules, 
fullerenes, etc. which in general agree with experimental data and DFT calculations. It 
can be used for simulations of, e. g., large organic molecules, clusters, and bulk carbon- 
based materials in the cases that ab initio calculations are limited by abilities of modern 
computers. 

A serious drawback of the model [32] is, however, unsatisfactory description of C-C 
bonding in systems with low coordination numbers of carbon atoms. For example, 
the calculated binding energy of the dimer C2 is almost a factor of two greater than 
its experimental value, while the bond length C-C is much shorter. For the chains C3 



and C4, the discrepancy between the theory and experiment is not so large but still 
rather appreciable (the binding energies are greater than the experimental values by ~ 
10% and « 20% respectively). Thus, the model overestimates the bond strength for 
carbon atoms with Kc = 1 and/or 2 in the absence of adjacent hydrogen and/or oxygen 
atoms. As a consequence, this model gives qualitatively incorrect results in the cases 
when accurate values of binding energies of small carbon clusters are needed (e. g., the 
fragmentation of fullerene Cgo through the loss of a C2 dimer [211 ) or for the systems with 
large relative number of carbon atoms having K(,= l and 2. For example, the 100-atomic 
graphene fragment appears to be unstable and loses its overall hexagonal structure upon 
separation of boundary atoms during relaxation (this does not occur if the model [H] 
developed for carbon systems is used). A possible reason for the mentioned shortcomings 
of the model [22] is that the parameters of the tight-binding Hamiltonian were fitted to 
DFT calculations and not to the experiment. 

The present paper is aimed at the search for such a set of parameters of the nonorthog- 
onal tight-binding model that does not result in contradiction with experimental data for 
the C2 dimer and a number of other carbon systems, while the energetical and structural 
characteristics of various CnHm molecules and clusters agree with experiment no worse or 
even better than in Ref. [22] • We fit the model parameters based on the criterion of the 
best correspondence between the calculated and experimental (not DFT-derived) values 
of bond lengths and binding energies of several selected small CnHni molecules. The 
resulting tight-binding potential appears to work well for other, relatively large CnHm 
molecules and clusters as well. It also correctly describes the crystalline carbon struc- 
tures. In this work, we restrict ourselves to hydrocarbons. Generalization to the C-H-O 
systems will be done later. 

The paper is organized as follows. In Section 2, we recall the nonorthogonal tight- 
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binding model for hydrocarbons, describe the alrorithm used to find its parameters, and 
present the results obtained. In Section 3, we compare the calculated bond lengths and 
binding energies of various CnHm molecules and carbon clusters with the corresponding 
experimental values and the results of the work . In Section 4 we consider bulk carbon 
structures, diamond and graphene, as well as hydrogen interstitial defects in diamond. 
Section 5 concludes the paper. 

2. Parameters of nonorthogonal tight-binding model for hydrocarbons 

In the tight-binding model, the total energy E for a given set of atomic positions {R^} 

is 

E = Eei + Ej-ep, (1) 

where 

^el= E (2) 

n,CT(occ) 

is the quantum-mechanical electronic ( "band" ) component of E which is the sum of one- 
electron energies for the occupied states (in the absence of magnetic field, e„ does not 
depend on the spin projection cr = j or j), 

^rep = EE0(i^^.) (3) 

i j>i 

is the classical component of E", being equal to the sum of pairwise ionic repulsive po- 
tentials = 0(|Ri — Rjl). 

The energy spectrum {£„} is found from the solution of the stationary Schrodinger 
equation 

H^nir) = £„*n(r) (4) 

by the expansion of eigenfunctions 

*n(r) = J:Cra^^a{r) (5) 

i,a 
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in terms of the nonorthogonal atomic orbitals {(^^^(r)}, where i is the atomic number, a 
labels the type of atomic orbital (IS* orbitals of hydrogen atoms and 2S, 2Px, 2Py, 2P^ 
orbitals of carbon atoms are taken into account). Upon substitution of Eq. (5) into Eq. 
(4), multiplication by ^jjsi^) from the left side and integration over r, equation for the 
one-electron self-energies takes the form 

UHi^ - enSi^)Cl = 0, (6) 

i,a 

where 

Hi^ = J drip*^{r)HipUr) (7) 
are the matrix elements of the Hamiltonian, 

Si^ = I drip]p{r)ipi^{r) (8) 

are overlap integrals for atomic orbitals. The one-electron Schrodinger equation (4) is 
thus reduced to the generalized eigenvalue problem (6) which is solved numerically. The 
number of equations in the system (6) equals to the total number of atomic orbitals 
involved, i. e., A^h + 4A^c, where A^h and A^c is the number of hydrogen and carbon 
atoms respectively. Molecular orbitals (5) are occupied (according to the Fermi-Dirac 
distribution function and Pauli principle) by electrons, the number of which is Ah + 4A^c 
as well. Note that in the orthogonal tight-binding models, matrix is diagonal in (z, j) 
and (q;,/?). 

If the nonorthogonal tight-binding model is used in molecular dynamics simulations, 
the force = —dE/dHk acting on the k-th atom is found after calculation of eigenen- 
ergies Sn and eigenvectors C-^. As follows from Eq. (6), 

t,aj,l3 



Explicit expressions for Sjj^ and (see below) greatly simplify simulations of relatively 
large systems. 
For Hfj^ the following parametrization is used 



Hl^ = -K,,Sr:,{H^ + Hp)., 



1 



(10) 



where 



1, ^ = j 

Kf- exp [-5^j [R 



V 



(note that H}^ = at a = [3 and H}'^ = at a 7^ since S}^^ = ^a/?)- There are 
three different parameters i^Q,, according to the number of different atomic orbitals a = 
IS, 2S, 2P, and three different parameters Kfp 6ij, R^j each, according to the number of 
pairs of atoms of different sorts ij = HH, HC, CC (note that in Ref. [22], the coefficients 
5ij were taken to be the same for all pairs of atoms). The values of Sj^ are obtained 



ri(3 



if3 



analytically [251 ES] using the expressions for Slater orbitals overlaps: 

exp (-6^0, 
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(13) 
(14) 



where j = x,y, z. 
Pair potentials in Eq. (3) for E'j.gp are taken in the form 

0(i?,,) = 0?,exp {-^,,{R,,-R% 



(15) 



(there are three parameters 0^- and i^ij each). So, the total number of fitting parameters 



(i^Q, A'J-, b^j, R^^^, 0^^-, (3ij, ^is, 6p) is 21. We remain unchanged (see Ref. [22j) 
the parameters His = -10.70 eV, K^^ = 1.68 eV, (5hh = 0.13 A^^, R% = 0.75 A, 
^HH = 0.78 eV, /3hh = 6.84 A~\ and ^is = 2.456644 A^^ describing the H-H interaction. 



As mentioned above, our main purpose was to considerably improve the correspondence 
between theory and experiment for small carbon clusters and several other purely carbon 
systems. We were unable to do it through changes in the parameters of C-C interactions 
only, because of drastic decrease in accuracy of the results obtained for various hydro- 
carbon molecules CnHm. This is why we changed the parameters of C-H interactions 
along with the parameters of C-C interactions. Search for new values of those parame- 
ters was based on the criterion of the best correspondence between the calculated and 
experimental interatomic distances, binding energies, and some vibration frequencies of 
the following clusters, molecules, and radicals: C2, C3, C4, CH, CH2, CH3, CH4, C2H2, 

Since the number of fitting parameters exceeded the total number of model parameters, 
the latter were derived by minimization of the function ^{H2S, H2P, Kqq, i^cH? •••) which 
we chose as a quadratic form of the differences between the corresponding calculated 
and experimental values. Usual minimization techniques such as the gradient or Newton 
methods can hardly be applied to the function $ because 1) this function has many 
local minima and 2) cusps or even discontinuities in the dependences of the physical 
quantities on the model parameters are possible, thus making the use of those methods 
questionable since they imply that the function $ is sufficiently smooth. For this reason, 
in order to find the global minimum of the function $ we used the simplest version of 
the Monte-Carlo method. Our strategy was to explore the neighbourhood of the last 
"best" (i. e., having the minimum value of $) point in the parameter space among the 
points considered up to a given iteration step. The coordinates of each new point were 
determined using the generator of pseudo-random numbers. The size of the region for 
the search of new points was periodically changed, which is necessary both for refining 
the coordinates of the $ minimum found at the preceding step and for finding other, 
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probably more deep minima. 

Finally we obtained the following values of the model parameters: H2S = —16.157972 
eV, H2P = -10.078261 eV, K^c = 2.060290 eV, K^^ = 1.763801 eV, ^cc = 0.164262 
A~\ Sen = 0.014350 A-\ R^cc = 1-582565 A, R^cu = 1-045120 A, c^^^c = 0.943505 eV, 
00^ = 0.561102 eV, Pec = 4.912617 A^^, Pen = 9.433587 A^^, ^2S = 2.991164 A'^, 
^2P = 3.857861 A~^ They do not differ much from those given in Ref. [32j, with the 
exception of 6cr which is a factor of ^ 9 smaller that in Ref. [321 • We draw attention to 
the fact that ^25* < ^2P, while ^25 > ^2P in Ref. [22 • 



3. Binding energies and structures of C^H^ molecules 

The binding energies Ei, of clusters and molecules CnHm were determined as 

Eb (C,H,,) = nE{C) + mEiR) - E(C,H^) , (16) 

where E (Cr^H„^) is the total energy of the system, E (C) and E (H) are the energies of 
isolated carbon and hydrogen atoms, respectively. In the tight-binding model, they are 
E{C) = 2£'25 + 2E'2P and ECR) = Eis- The specific (per atom) value of E^ is found from 
Eq. (16) through deviding by the total number of atoms in the system, (n + m). The 
results obtained for several clusters and molecules are hsted in Table 1 along with the 
corresponding experimental values and the values calculated using the set of the model 
parameters from Ref. [22] • Note that theoretical values of Ei, are given without account 
for the zero-point energy which IS, as a rule, about ~ 0.1 eV/atom and results in a small 
decrease of E^. 

From Table 1 one can see that the binding energies for the dimer C2, the trimer C3, and 
the chains C4 and C5 are much closer to the experimental values than those calculated 
with the model parameters from Ref. [22] (for the dimer C2 there is a qualitative change 
in the value of E^ almost by a factor of two as compared with Ref. [|32j). For the 
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majority of small CnHm molecules and CgHg cluster our values of Eb are also closer to 
experimental ones than those calculated with the parameters of Ref. [32], see Table 
1. This is true for relatively large molecules CnHm as well. For example, we obtained 
Eb = 5.09; 4.31; 5.03 eV for CioHg (naphthalene), CioHie (adamantane), and C12H10 
(acenaphthene) respectively, while the experimental and calculated with the parameters 
from Ref. |32j values are, respectively, Ej, = 5.07; 4.35; 5.03 eV |S7| and Et = 5.17; 4.40; 
5.09 eV. 

Table 1 shows also theoretical and experimental values of interatomic distances C-C 
and C-H. One can see that for the dimer C2 there is much better correspondence with 
the experiment as compared with Ref. [32]. For the molecules and large clusters CnHm, 
the deviations from the experimental values, as a rule, do not exceed several hundredth 
of A, such an accuracy being sufficient in simulations of a broad range of hydrocarbons. 
However, it should be pointed out that the C-H bond lengths calculated in Ref. [22 



are 



often somewhat closer to the experimental ones. As mentioned above, this is at a price of 
unsatisfactory description of C-C interactions in some carbon systems, while our model 
adequately describes, e.g., the two-dimensional graphene fragments which preserve their 
overall hexagonal structure upon relaxation, being just shghtly distorted. 

We have calculated the binding energies and bond lengths in several fuUerenes. For the 
smallest possible fullerene C20 [22], we obtained Ei, = 6.31 eV/atom, in agreement with 



the values of E^ = 6.08 eV/atom [IB] and E^ = 6.36 eV/atom [40] found, respectively, 
within the orthogonal tight-binding model IH] and by DFT with gradient corrections (we 
are not aware of experimental Ei, value for the fullerene C2o)- The calculated minimal 
and maximal C-C bond lengths, /^m = 1.44 A and Imax = 1-52 A coincide with those 
obtained within the orthogonal tight-binding model [[H] and agree well with recent ab 



initio calculations for various exchange-correlation potentials {Imin = 1-40 — 1.43 A, 
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Imax = 1-51 — 1.52 A). Note that in the model [32j, the values of = 5.89 eV/atom, 
Imin = 1-45 A, and Imax = 1-61 A differ from the results of Refs. ^B, HDl Ej much 
greater. 

As for the different isomers of C20, we have found that the bowl isomer is energetically 
more favourable than the cage, its binding energy being Ei, = 6.55 eV/atom, while the 
ring isomer with Ei, = 6.81 eV/atom is more stable than the bowl one. In the absence 
of experimental information on the relative stability of C20 isomers, various theoretical 
approaches give conflicting results. However several authors arrived at the same sequence 
of C20 isomers as described above, see references in Ref. 

Next we calculated the coagulation energy AE = 2£'[C2o] — ^[(^20)2] of two fullerenes 
C20 in the open-[2+2] isomer of a cluster molecule (€20)2 [42]. We found AE" = 4.9 eV. 
This value coincides with that obtained in the orthogonal tight-binding model jH, 



and does not differ much from the value of AE" = 6.3 eV obtained in the DFT [ i42j . 
The intercluster bond length / = 1.37 A agrees well with both orthogonal tight-binding 
model HH (/ = 1.35 A [22j) and DFT (/ = 1.34 A [42]). 

For fullerene Ceo 113] , the binding energy and the bond lengths are E^ = 7.01 eV/atom 
and / = 1.41, 1.48 A, respectively, in good agreement with the experimental values of 
Eb = 6.93 eV/atom and / = 1.40, 1.46 A [37]. The model [32] gives E^ = 6.71 eV/atom 
and / = 1.42, 1.53 A, in a poorer correspondence with the experiment. Meanwhile, in 
the model [32], the HOMO-LUMO gap A = 1.91 eV is closer to the experimental value 
A = 1.6 — 1.8 eV [H] than our result A = 1.15 eV. So, the suggested set of parameters of 



the nonorthogonal tight-binding model is more suited to calculations of binding energies 
and bond lengths than electron characteristics. 

For fullerene C70 we obtained E^ = 7.04 eV/atom and / = 1.41 — 1.49 A, again in much 
better agreement with the experimental data E^ = 6.97 eV/atom [15] and / = 1.37— 1.48 
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A [IB] than the results of the model [32j = 6.73 eV/atom and / = 1.41 - 1.55 A (we 
are not aware of reliable experimental data on the heats of formation and bond lengths 
in fullerenes other than Cqq and C70). 



4. Bulk crystalline forms of carbon 

We have also calculated the binding energies and bond lengths for two crystalline forms 
of carbon, diamond and graphene. We made use of periodic boundary conditions and 
performed scaling with respect to dimensions of a supercell. For the binding energy and 
bond length in diamond we obtained = 7.36 eV/atom and / = 1.54 A, in excellent 
agreement with experimental values = 7.35 eV/atom and / = 1.54 A. On the other 
hand, the model [i32j gives E^ = 6.58 eV/atom and / = 1.61 A. These results differ 
substantially from experimental ones. 

For graphene our tight-binding potential results in Ef, = 7.36 eV/atom and / = 1.45 A. 
Since the experimental value of E^ in graphite is 7.37 eV/atom and the weak interlayer 
coupling is abouth 0.02 eV/atom, for graphene one has E^ = 7.35 eV/atom, in close 
correspondence with our value of Ei,. Although the calculated bond length is somewhat 
large than the experimental one, / = 1.42 A, we note that the model ||32] leads to much 
stronger deviations from the experiment in both £"5 = 7.06 eV/atom and / = 1.48 A. It 
is particularly remarkable that a rather good description of energetics and structure of 
bulk phases of carbon was given within the model whose parameters had been derived 
based on a comparison with experimental data for small clusters and molecules. This 
points to a good transferability of our tight-binding model, at least in what concernes 
pure carbon systems. 

Now we turn to a much more subtle problem, the hydrogen interstitials in diamond 
[ ITT] . Unlike Si and Ge, there is no experimental data on isolated interstitial hydrogen in 
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diamond. Based on the experiments with muonium, it is generaUy beheved that there 
are at least two stable interstitial sites, tetrahedral (T) and bond-centred (BC), the latter 
being 1 - 2 eV lower in energy (see Fig. 1 in Ref, [HT]). 

We made calculations for 64- and 216-atom supercells. The results vary only slightly 
with the size of a supercell. We have found that the BC-interstitial corresponds to a 
minimum of the total energy as a function of atomic coordinates, i.e., all oscillation 
frequencies are real. The formation energy of this interstitial is Ej = 1.4 eV. Two C-H 
bond lengths are both 1.10 A, in accordance with calculations of other authors, 1.05 
- 1.17 A jl^. We have also identified several other stable interstitial sites, including 
the H-interstitial with the energy 1.3 eV higher than that of the BC-interstitial (to be 
compared with the values 1.5 - 1.9 eV obtained within DFT [H7]). 

Detailed description of our results on hydrogen interstitials in diamond goes beyond the 
scope of this paper. We note, however, that the T-interstitial appeared to be the saddle 
point lying 1.4 eV above the BC-interstitial. Upon relaxation, the hydrogen atom moved 



to one of two stable sites close to the C-site [17]. The reasons why the T-interstitial 
appeared to be not a local minimum of the total energy are not quite clear to us. The 
only essential difference between the multi-coordinated hydrogen atoms in H- and T-sites 
is that in the latter case the bonds C-H do not lie in the same plane. 



5. Conclusions 

The set of parameters for the Hamiltonian of the nonorthogonal tight-binding model 
derived in this work not only provides rather accurate description of C-H interactions 
in hydrocarbons C„H„i, but also allows one to simulate the structure and energetics of 
both small (C2, C3, C4) and relatively large (fuUerenes C20 and Ceo, graphene fragments) 
carbon clusters, as well as bulk crystalline forms of carbon (diamond, graphene). The 
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accuracy of the method is sufficient for many purposes. Analytical dependences of the 
overlap integrals on atomic coordinates greatly simplify the calculation of forces acting on 
the atoms. This makes possible the molecular dynamics simulations of relatively large 
systems for which ab initio calculations are problematic due to the limited computer 
power. Besides, the evolution of small clusters can be followed on a microsecond time 
scale jlSlllS] (the characteristic times are about 1 - 10 ps for ab initio molecular dinamics). 



On the other hand, this tight-binding potential fails to describe the tetrahedral hydrogen 
interstitial in diamond, although the bond-centred and hexagonal interstitials are well 
reproduced. Hence it may be not applicable to the systems with multi-coordinated 
hydrogen atoms where the C-H bonds do not lie in the same plane. There are, however, 
just few examples of such systems. Later we plan to generalize our approach to the 
C-H-0 systems. 
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Table 1. Binding energies Eb and bond lengths /cc, ^ch for some Cnilm molecules. 

* Experimental values of E'^feV/atom] are obtained from known heats of formation 
AfH° [kJ/mol] [37| making use of relations Eb = j^y{nE^ + mE^ - 1.0364 • lO^^ . A///°}, 
where Eg = 7.3768 eV, E^ = 2.375 eV [9]. 
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Formula 
and 
Name 


Binding energy eV/atom 


Bond lengths, A 


Potential 
from [32] 


Present 
work 


Experimental 
data* [SZl EBj 


Potential 
from [|32j] 


Present 
work 


Experimental 
data pSj 


C2 
Carbon 
dimer 


5.78 


3.15 


3.12 


1.157 


1.230 


1.243 


C3 
L/arbon 
trimer 


5.10 


4.72 


4.54 


1.300 


1.301 


1.277 


C4 
Carbon 
chain 


5.86 


5.09 


4.88 


1.187 
1.513 


1.296 
1.354 


— 


C5 
Carbon 
chain 


6.07 


5.68 


5.35 


1.263 
1.338 


1.273 
1.348 


— 


CH 

Methyhdyne 


1.90 


1.87 


1.80 


1.089 


1.081 


1.120 


CH2 
Methylene 


2.80 


2.75 


2.69 


1.079 


1.080 


1.085 


CH4 
Methane 


3.66 


3.40 


3.51 


1.089 


1.100 


1.094 


C2H2 
Acetylene 


4.87 


4.54 


4.29 


C-C 1.199 
C-H 1.064 


C-C 1.226 
C-H 1.079 


C-C 1.203 
C-H 1.063 


C2H4 
rjthylene 


4.19 


3.96 


3.94 


C-C 1.315 
C-H 1.091 


C-C 1.327 
C-H 1.097 


C-C 1.339 
C-H 1.086 


C3H4 
AUene 


4.56 


4.32 


4.22 


C-C 1.310 
C-H 1.096 


C-C 1.323 
C-H 1.100 


C-C 1.308 
C-H 1.087 


CqRq 

Benzene 


4.95 


4.82 


4.79 


C-C 1.422 
C-H 1.092 


C-C 1.407 
C-H 1.095 


C-C 1.397 
C-H 1.084 


CgHg 
Cubane 


4.37 


4.42 


4.47 


C-C 1.645 
C-H 1.075 


C-C 1.570 
C-H 1.082 


C-C 1.571 
C-H 1.097 
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